Shear Zones in granular materials: Optimization in a self-organized random potential 
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We introduce a model to describe the wide shear zones observed in modified Couette cell experi- 
ments with granular material. The model is a generalization of the recently proposed approach based 
on a variational principle. The instantaneous shear band is identified with the surface that minimizes 
the dissipation in a random potential that is biased by the local velocity difference and pressure. 
The apparent shear zone is the ensemble average of the instantaneous shear bands. The numerical 
simulation of this model matches excellently with experiments and has measurable predictions. 
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Strain localization or shear band formation in granu- 
lar materials has been studied for many years Q, Q due 
to its importance in engineering and geoscience. Shear 
bands are narrow regions separating almost solid blocks 
moving with different velocities. Recently Fenistein and 
coworkers observed 0, Q wide shear zones in a modified 
cylindrical Couette cell. Instead of letting an inner cylin- 
der rotate with respect to an outer one, their cell has no 
inner wall, but exerts the shear deformation via the bot- 
tom, which is split into a rotating inner disk of radius R s 
and a fixed outer annulus (see Fig. The shear zone 
is pinned at the bottom split and evolves independently 
of walls. This approach has attracted considerable in- 
terest 0, 0, because it provides new insight into the 
fundamental problem of shear band formation. 

A theory based on the principle of minimum dissipa- 
tion rate was proposed soon after the first publications. 
With the assumption of negligible width of the shear zone 
a model with no fitting parameter resulted which was 
used to describe the position of the shear zones 5] . This 
model proved to be efficient and even delivered predic- 
tions about a new type of closed shear zones, which have 
been found later both in experiments and in computer 
simulations [f| Q. However, the model has to be gen- 
eralized in order to describe the interesting phenomena 
related to the width of the shear zone, and this is the aim 
of the present Letter. 

Let us first briefly summarize the experimental find- 
ings: The shear zone starts from the bottom split, and 
for small and moderate filling H < 0.7 R s it ends up on 
the surface 0. If the filling height is further increased, 
the shear zone is buried in the material and takes the 
shape of a cupola 0, 0, 0] ■ The surface position of the 
shear zone for small filling can be very well described by 
a universal empirical curve. The width of the shear zones 
on the surface increases as a power law with the filling 
height with an exponent ~ 2/3. For the shape and width 
of closed shear zones, however, there are only qualitative 
experimental results so far. 



Our model combines the ideas of two main sources. 
The first one is our former theoretical analysis of shear 
zones based on the principle of minimum dissipation 
That model (called "Unger-model" for short) assumes 
that the shear zone is infinitely narrow, separating the 
standing and moving parts of the material. The local 
dissipation rate is proportional to the velocity difference 
across the shear band and to the hydrostatic pressure. 
In this simple model the position of the shear band is 
identified with the shape that corresponds to the global 
minimum dissipation rate. By definition, this model was 
only able to describe the position of the shear band. Due 
to the cylindrical symmetry the problem was traced back 
to finding the minimum path in a smooth 2D potential. 
In spite of its simplicity, the model gave surprisingly ac- 
curate results, moreover, it predicted the transition from 
the open (Fig. [I]b) to the closed shapes (Fig. |T]c). 

Another optimization problem was introduced earlier 
in the context of shear band localization during compact- 
ification in sheared loose granular matter |8(. In that 
model an instantaneous shear bands were identified by 
the global minimum of a random potential representing 
the local inhogeneities of the granular material. These 
inhogeneities are not a priori frozen in but change due 
to the relative displacement in the shear band. Thus the 
global minimum may be at a different place in the next 
moment giving rise to an ansamble of instantaneous shear 
bands which can be identified to the visible macroscopic 
displacement field. These two models are here combined, 
and the fluctuations induced by the random potential 
are used to address the problem of the width of the shear 
bands. 

We generalize the Unger-model in the following wa~ 
Instead of using a smoothly varying potential like in 
we use a random one which is suited better to the disor- 
dered nature of granular media. Like in ref. the in- 
stantaneous displacement is represented by a single local- 
ized shear band determined by optimization on a random 
field: The shear band is the path which minimizes the dis- 
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FIG. 1: a) The experimental arrangement and caluclated 
shear bands. The dotted circle is the bottom split, the defini- 
tion of the most important quantities are noted on the figure. 
The gray scaling on all plots is proportional to the logarithm 
of the occurrence probability at the given point. R a — 150 
for a-c) and H = 115, 100, 150 for plots a, b, c) respectively. 



sipation rate and fits the boundary conditions. (We shall 
again assume that, due to the symmetry of the problem, 
the determination of the band reduces to that of a line.) 
The shear is known to change the local structure of the 
material which we take into account by changing the ran- 
domness in the neighborhood of the actual shear band. 
A new optimal shear band is then searched for. For each 
random realization of the material the minimization de- 
termines a single narrow band. The shear zone itself is 
represented as an ensemble of narrow bands, i.e. the flow 
velocity of the material can be obtained as an ensemble 
average over the realizations. Our approach is related 
to the first passage percolation problem |j, also known 
as polymer in a random medium 10], with the difference 
that in our case the randomness organizes itself dynam- 
ically. 

The model is defined on a regular square lattice which 
is applied as a coarse grained representation of the mate- 
rial in a radial cut from the center to the outer wall. In 
the radial cut the shear band is represented by a contin- 
uous path that starts from the split point of the bottom 
and reaches either the surface or the axis of the sample. 
We allow nearest and next nearest neighbour connections 
with lengths A£ equal to 1 and v2 respectively. Such a 
path V stands for a possible sliding surface with cylin- 
drical symmetry. 

The energy dissipation rate associated with path V is 
in this geometry proportional to the torque due to the 
local friction forces. These are modelled by a random 
strength parameter u(r, h) assigned to each site of the 
lattice, where r and h are the radial and height coordi- 
nates. Each random variable u is generated uniformly 
in [au max , u max ], < a < 1. As in Ref.[f| we set 
Umax = f 2 (H — h) , which up to constant factors is hy- 
drostatic pressure (H — h) times cylindrical circumfer- 



ence (2?rr) times lever (r). The torque is then given by 
integrating the local shear resistance of the material over 
the path V: S = J2ie-p Ui ^ i - The actual instantaneous 
shear band is the directed path that can be activated 
by the smallest torque, i.e. the one which minimizes S. 

Once the minimal path is found we refresh the strength 
parameters u randomly along it and in its vicinity (near- 
est neighbour sites). By the successive application of this 
procedure an ensemble of shear bands is collected which 
provides the velocity field of the shear flow. One instan- 
taneous shear band separates the sample into two parts 
in such a way that each site in the inner part rotates by 
the driving angular velocity Q, of the bottom disk while 
the sites in the outer part have angular velocity zero. 
Taking an average over the ensemble of shear bands one 
arrives at a field of angular velocity that can be compared 
directly to that observed in experiments. 

An important difference between the model of [f| and 
the present one is that the randomness is refreshed not 
only on the current minimal path but also in its neigh- 
borhood. If only the values in the path were changed 
no steady state would be reached. The average value of 
the randomness u would increase continuously, however 
extremely slowly. This is not the case here where the 
steady state is reached fast as in experiments. This vari- 
ant of the model is closely related to the self-organized 
criticality model of Bak and Sneppen 0] • 

The model has three free parameters. In dimension- 
less quantities they are: R s /H, R s /a, a. In words: the 
aspect ratio of the sample, the split radius in units of 
the lattice constant and a number, which controls the 
effect of disorder in the model. Henceforward in this pa- 
per we use a— I, i.e. all lengths are measured in units 
of a. The parameter a mimics strength fluctuations due 
to individual properties (shape, friction, etc.) and coop- 
erative effects (e.g., density fluctuations). A smaller a 
means stronger disorder, while a = 1 would be a lattice 
version of the deterministic Unger-model|l3j. We present 
data for several different values a = 0, 0.43, 0.5, 0.6 and 
many different split radii in the range of R s = 15 . . . 600. 
As the lattice unit must be larger than the lower length 
cutoff provided by the particle diameter, larger values of 
R s correspond to smaller grain size. 

In Fig. ^the probability distribution of instantaneous 
shear band positions is plotted. We get very similar pat- 
terns as the ones obtained in experiments and molecular 
dynamics simulations 0, Q with both open and closed 
shear zones. We calculate the angular velocity at any 
point of the sample and compare it to the experiments. 

If the shear zones are far from the system boundaries 
the error function is a very good fit for the angular veloc- 
ity in agreement with the experiments It gives both 
position and width of the zones. 

Most of the experimental data concern the surface po- 
sition of the shear zones (R c )- We compare first this 
property on Fig. [21 The analytical result of p| overesti- 
mated the experiments for H/R s > 0.25. Increasing ran- 
domness decreases the apparent shear zone radius with 
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FIG. 2: The surface positions of the shear zones. The solid 
line is from the dashed line is the experimental 3] curve. 
Symbols were obtained for systems with R s — 90 and a = 
0, 0.43, 0.5, 0.6 for +, x, o, □, respectively. The inset shows 
the same data on log-log plot. 




FIG. 4: The scaled center position (curves a) and width 
(curves b) of the closed shear zones. The solid line is the cal- 
culated position from|f|. Datapoints (+, x, o, □) connected 
with dashed lines have i? s =90, 150, 300, 600 and a=0.43. 




FIG. 3: The surface width of the shear zones. The solid 
line is {H / R s ) 2 ^ ■ The lower curves are for a = 0.43 with 
R s = 90, 150, 300, 600 with symbols +, x, *, □, respectively. 
The upper curves are for a = with R s = 150, 300 with 
symbols o, A respectively. The breakdown of scaling close to 
the transition is shown in the inset for a — 0.43. 



the best matching at about a ~ 0.5. System size does 
not influence the curves for R s ^> 10. 

From the experiments, Fenistein et al. 3] suggested 
a power law dependence of the surface position on the 
height, namely 1 - R c /R s = (H/R s ) 2 5 . Thus on the 
inset of Fig. [5] we tested it on a log-log plot. Both the 
Unger-model calculation and our numerical data devi- 
ate slightly from the simple power law function for small 
H/R s . This deviation is too small to be seen on a nor- 
mal plot but could be tested on the experiments if high 
enough precision can be attained [l4l ]. 

In the experiments the width of the shear zones on 
the surface was found to be a power law of H with an 
exponent of 2/3 Directed polymers have the same 
roughening exponent |10|. As shown on Fig. [3] we ob- 
tain an exponent very close to this value. The curves for 
different R s (but the same a) can be scaled together by 

plotting W/R^ 3 versus H/ R s as in the experiments . 

This power law increase of W with H must stop, when 
the width of the shear zone reaches R c , i.e. the available 
distance between the container axis and the average shear 
band position at the surface. For larger H this finite size 



effect implies that W w R c cx R s (for given H/R s ), which 
explains the sudden increase and the loss of data collapse 
of the scaled curves in Fig. |3|at about H/R s ~ 0.7. This 
feature can also be observed on the experimental results 
0. It cannot be interpreted as a sign of an increasing 
characteristic length scale at the transition from open to 
closed shear zones. 

We find that with decreasing a (i.e. increasing dis- 
order) the width increases. This agrees with the exper- 
iments, where more irregular particles produced wider 
shear zones Q • If a is varied, not only the surface width 
but also its position changes slightly. This seems to be 
a secondary effect not yet observed experimentally. Our 
prediction is that the scaled surface position R c /R s of 
the shear zones is slightly larger for smoother particles 
than for rough ones. 

Due to the occurence of closed shear zones the local an- 
gular velocity on the symmetry axis of the container de- 
pends on h, being equal to fl at the bottom and decreas- 
ing monotoneously towards the surface. By fitting this 
dependence again by an error function, we determine the 
position htop and vertical width W c of the closed shear 
zones. This works as well as on the surface provided the 
shear zone is not too close to the boundaries. This fitting 
procedure has a much broader range of applicability than 
the one with half a Gaussian which works well only for 
small systems R s < 30 with moderate H < R s |7J. 

Some of our measured datasets are plotted on Fig. 0] 
The closed shear zones become flatter with increasing H. 
The position scales with R s and follows the curve calcu- 
lated in p| for values of H/R s between 0.8 and 2. The 
deviations for larger H/R s can be understood by noting 
that the height of the shear zone decreases faster than the 
width with increasing H. Above a certain filling height 
the shear zones touch the bottom of the container. This 
raises the apparent position of the shear zones compared 
to the noisless system of Q- 

The vertical width W c of the closed shear zones scales 

2/3 

with R s , as expected for directed polymers with a 
length of the order R s . Close to the transition, for H/R s 
between 0.7 and 0.8, both W c and h top deviate from the 
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FIG. 5: The order parameter m for q = 0. System size R 3 — 
15, 35, 75, 90, 150, 300, 600 from right to left respectively 
(+, x, *, □, o, A, +). Experimental data from 6] are 
shown with filled symbols of 0, v, □ for R a = 45, 65, 95cm 
respectively. The inset shows the width of the tranisition 
versus R s for numerics + and experiments 0. The dashed 
line has the slope of 0.5. 

expected behaviour. Again this can be understood as a 
boundary effect: In the absence of fluctuations, h top is 
closer to the surface than the width W c permits. 

The next question we focus on is the phase transition. 
There seems to be a discrepancy between the theoretical 
model which predicts a first order and the experiments 
0,0 which claim to see a continuous transition. The lat- 
ter two experimental papers introduce different empiri- 
cal fittings of the shear zone profiles and show the tran- 
sition of the fitting parameters. Both papers estimate 
that the transition occurs at lower H/Rs (H/Rs — 0.6 
0, H/R s ~ 0.65 H) than in the theory of 5]. The 
drawback of both approaches is that only one side of the 
transition can be studied. 

We prefer the classical approach of the order parame- 
ter of the transition. A good candidate seems to be the 
normalized angular velocity of the surface at the center 
m = uj(r—0, h=H)/fl. If the system has only open shear 
bands: m = 1, if only closed ones: m = 0. If both types 
are present, m can take any value between and 1. 



Fig. [S] shows the change of m with H/R s for different 
R s and a = 0. The transition gets sharper as the system 
size increases: In the thermodynamic limit R s — > oo, 

— 1/2 

its width seems to vanish like R s (see inset of Fig. 
|SJ. Then the order parameter jumps at a value of H/R s , 
which wc could estimate from a finite size scaling analysis 
as (H/R s ) c ~ 0.735. Remarkably, this is very close to the 
higher limit of the hysteresis calculated in [5| . 

The angular velocity of the surface at the center is 
available in experiments so that the test of the order 
parameter is straightforward. The plot of m for three 
systems can be found in Ref. [6j. The experimental re- 
sults look quite similar but with a little shift in H/R s . 
The sharpening of the transition in experiments cannot 
be tested due to the limited range of R s . 

In conclusion we have shown that the results obtained 
from numerical simulations of our model can be directly 
compared to the experimental ones. Excellent agreement 
can be obtained for the already measured quantities such 
as surface position, width and angular velocity at the 
center of the surface. The comparison with experiments 
shows that our lattice constant can be indcntificd roughly 
with one particle diameter. This also draws attention to 
the fact that the experimental systems R s = 15... 95 
are far from the thermodynamic limit especially if one 
studies the order of the transition. More attention should 
be given to this fact. 

The variational principle combined with the self- 
organized random potential appears to be an efficient 
tool to study shear zones maybe also in other geome- 
tries. The simplicity of the model with just a single con- 
trol parameter for the local strength fluctuations makes 
it robust, and applications to different geometries seem 
to be straightforward. 
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